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ABSTRACT 



A region of a star that is stable to convection according to the Ledoux cri- 
terion may nevertheless undergo additional mixing if the mean molecular weight 
increases with radius. This process is called fingering (thermohaline) convection 
and may account for some of the unexplained mixing in stars such as those that 
have been polluted by planetary infall and those burning ^He. We propose a new 
model for mixing by fingering convection in the parameter regime relevant for 
stellar (and planetary) interiors. Our theory is based on physical principles and 
supported by three-dimensional direct numerical simulations. We also discuss the 
possibility of formation of thermocompositional staircases in fingering regions, 
and of their role in enhancing mixing. Finally, we provide a simple algorithm 
to implement this theory in one-dimensional stellar codes, such as KEPLER and 
MESA. 

Subject headings: Convection, Diffusion, Hydrodynamics, Instabilities, Planet-star 
Interactions, Stars: Evolution 
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1. Introduction 

The phenomenon of fingering convection (otherwise known as thermohaline convection) 
occurs in regions of stellar and planetary interiors where the mean molecular weight, fi, 
increases upwards but which are nevertheless stable to the Ledoux criterion. If not for the 
effects of diffusion, such a system would be stable to small perturbations. In the presence 
of diffusion, however, instability can occur. Indeed, when a high fi, high entropy parcel is 
displaced downward, heat rapidly leaks via thermal diffusion. Since compositional diffusion 
is much weaker, the parcel remains denser than its surroundings and proceeds to sink. The 
instability takes the form of finger-like structures and can significantly increase the flux of 
heat and composition above that of molecular or radiative diffusion. Fingering convection 
appears in several key problems in astrophysics in which an inverse /i-gradient is observed, 
the most notable of which are planet infall and ^He fusion. 

Take for instance the problem of planet infall. Stars with detected planets have 



higher metallicities than those without (e.g. Santos et al. 2001; Fischer & Valenti 2005). 



The relevant question is then whether this observation results from a superior ability of 
high-metallicity gas to form planets or from pollution by planet infall. In the latter scenario, 
it is thought that stars accrete planets onto their outer convective zones, raising the 



metallicity at the surface. Vauclair (2004) was the flrst to argue that flngering convection 
may play an important role in this problem. Because the mean molecular weight in the 
convective zone rises upon absorbing a planet, the radiative region just beneath experiences 
an inverse /x-gradient and becomes unstable to the flngering instability. However, the 
time scale of the mixing by flngering convection was until recently unknown but is needed 
to determine how long the post-infall excess metallicity of the convective zone remains 
observable. Thus, this time scale determines whether the planet-metallicity connection is 
an effect of planet formation or infall. 
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During the process of ^He fusion, in which ^He + ^He — > ^He + 2^H, the total number 
of particles increases, so the mean molecular weight decreases ( Ulrich|[l972 ). Because fusion 
occurs more frequently at higher temperatures and densities, this reaction preferentially 
decreases the molecular weight deeper in the star, and builds up an inverse /i-gradient. 
This is thought to occur in red giant branch stars. Indeed, these stars are observed to 
undergo a process known as dredge-up in which particular isotopes, such as ^He, produced 



Gilroy 


1989 


Charbonnel 


1994 



Simulations by Charbonnel ( 1994 ) have produced surface abundances that agree well with 
observations for the so-called "first dredge-up;" however, she noted that some additional 
mixing in low-mass red giants is needed to account for the observed carbon isotope ratios. 



Charbonnel & Zahn (2007) later proposed that this additional mixing could be due to the 
fingering instability after ^He burning but were hampered by the lack of constraints on the 
mixing efficiency of fingering convection. Further progress on both subjects — planetary 
infall and ^He fusion — requires the development of improved mixing theories. 

Although the applications in which we are interested are astrophysical, much of 
the formalism of fingering convection was first developed to address the phenomenon of 
salt-fingers in the ocean. There, heat and salt play the role of the competing stratification. 



where salt diffuses about 100 times more slowly than heat. Stern (1960) first proposed the 
theory of a "salt-fountain" to explain the persistent, small-scale motions observed in warm, 
salty water lying above cool, fresh water even though density decreases upwards. For a 



recent review of the field, see Kunze (2003). By contrast with astrophysical fingering, it is 



possible to run laboratory experiments of fingering convection in salt water and measure 



its turbulent transport efficiency (e.g. Turner 1967 Stern & Turner 1969 Schmitt 1979). 
Unfortunately, none of these laboratory experiments are applicable to the parameter regime 
relevant to astrophysics. Because of this, theoretical models of fingering convection in 
astrophysics have remained, until recently, mostly phenomenological. 
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Several prescriptions have been suggested over the past four decades. Ulrich (1972) and 



Kippenhahn et al. (1980) attempted to constrain the dimensions of fingers using stabihty 
arguments and determined the resulting thermal and compositional mixing timescales using 
dimensional analysis. Both these prescriptions have a free parameter: in the former, the 
free parameter is the "effective inertia of the flow;" whereas in the latter, it is the aspect 



ratio of the fingers. Later, Schmitt ( 1983 ) used linear theory to predict the ratio of the heat 
to the compositional turbulent fluxes in astrophysical fingering but could not determine the 
absolute fluxes directly from this method. 

Only in the last few years have numerical simulations of fingering convection 
approaching the astrophysical parameter regime become more readily available (e.g. 



Traxler et al. 


2011a 


)• 


Denissenkov 


(2010 



the aspect ratio of fingers and used the previous literature and a dimensional argument 
to find a prescription of the mixing caused by fingering convection in red giant branch 
stars. He concluded that while this process could provide some of the required mixing 



discussed by Charbonnel &; Zahn (2007), it alone was insufficient to account for the 



observed abundances of low-mass red giant branch stars. Traxler et al. (2011a) presented 
3D numerical simulations of fingering convection and generated an empirical fit of their 
results to propose transport laws for fingering convection in astrophysics. Using their 



mixing model, Garaud (2011) studied the evolution of the surface metallicity of solar-type 
stars after planet infall and found that fingering convection would transport the material 
out of the convective zone too quickly to explain the planet-metallicity connection. This 
then suggests that planets can from more easily in high-metallicity proto-planetary disks. 

However, much work remains to be done. The failure of existing fingering convection 
models to explain red giant branch star abundances suggests that our understanding of 



the problem remains incomplete. Furthermore, as noted by Vauclair & Theado (2012), the 



- 6 - 



mixing prescription of Traxler et al. (2011a) does not fit their data well for systems which 



are only weakly stratified. Based on these considerations, our work here has several goals. 
First, we shall extend the available experimental datasets by running additional simulations 
closer to the true astrophysical regime. Second, we shall develop a theoretical — rather 
than empirical — prescription for the turbulent fluxes of heat and composition in fingering 
convection. And finally, we shall investigate other mechanisms of transport in fingering 
regions that may account for the additional mixing needed in the red giant branch case. 

Several such mechanisms have been proposed in the oceanographic literature. Not long 



after the initial theory of salt fingers was proposed. Stern (1969) realized that fingering 



convection had the peculiar property of driving coherent, large-scale gravity waves by an 



instability he termed the "collective-instability," which Stern &; Turner (1969) confirmed 



experimentally and which we address in more depth in § |5.2[ Another possible large-scale 
outcome of fingering convection in the ocean is the generation of thermohaline staircases 
(c.f. observations in the tropical Atlantic by Schmitt et al.||2005 ). These staircases are stacks 
of distinct, well-mixed, convective layers separated by thin fingering interfaces. Currently, 
the most promising explanation for staircase formation is the 7-instability, proposed by 



Radko (2003) and supported by numerical simulations (Stellmach et al. 2011). Both 



large-scale instabilities mix material much more efficiently than "homogeneous" fingering 



convection in the ocean (Schmitt et al. 2005). The 7-instability and collective-instability 



theories can be applied to the astrophysical case with minor modifications (for details, see 



Traxler et al. 2011a). However, because the viscous and compositional diffusion scales are 



minuscule in stars, it is currently difficult to perform full 3D numerical simulations in the 
astrophysical regime that can resolve both the fingers and these large-scale instabilities. 



From a theoretical point of view, Traxler et al. (2011a) predicted that layers could not form 



by the 7-instability in astrophysics, but we now investigate this more thoroughly and also 
pose the question of whether layers may form by the collective-instability instead. 
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In § [2| we describe our model. In § [3| we present the results of existing and new 



numerical simulations and compare them to the model described in Traxler et al. (2011a); 
we find that the latter does not fit our data as the system approaches overturning convection 



as suggested in Vauclair & Theado (2012). In § |4j, we then provide a new model that more 
precisely fits the simulations and that is constructed from physical principles rather than 
empirical fits. We discuss the conditions of layer formation in § [5] and conclude in § [6] by 
discussing this new model and its implications in stellar astrophysics. 



Model 



As discussed in Traxler et al. (2011a), the characteristic length scale of the fingering 
instability in astrophysical objects is much smaller than a pressure scale height, which in 
turn is generally small enough to ensure that the curvature of the star plays little role 
in the system dynamics. We therefore use a Cartesian grid (x, z), with gravity given 
by g = —gSz to model a small region of the star. In addition, the velocities within these 
fingers are much slower than the sound speed of the plasma, so we can treat the fiuid 



with the Boussinesq approximation (Spiegel & Veronis 1960). We assume that there is a 



background temperature profile, To{z), and a mean molecular weight profile, fio{z), which 
depend only on z. If the region considered is small enough, these background profiles 
can be approximated by linear functions with constant slopes, Tq^ and /ioz, respectively. 
We therefore assume that each profile is comprised of a linear background state and a 
perturbation, e.g. Ttot = zTq^ + T. 

The equation of state for the density perturbation, within the Boussinesq approximation, 

is 
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^ = -aT + /3/i, (1) 
Po 

where po is the mean density of the region, a = —p^^dp/dT, and (3 = p^^dp/dp. The 



remaining governing equations for a Boussinesq system are (Spiegel & Veronis 1960) 

V-u = 0, 



(2) 



— + u-Vu = ^+ g{aT - (3p)e, + z/V'u, (3) 

at pq 



dT 
~dt 



+ u-VT + w{To,- Tf) = ktV^T, and 



dp 

'di 



+ u-Vp + wpoz = K^V^/i, 



(4) 
(5) 



where u = {u, v, w) is the fluid velocity, p is the pressure, v is the kinematic viscosity, 
kt is the thermal diffusivity, /t^ is the compositional diffusivity, and T^'^ is the adiabatic 
temperature gradient, assumed constant within the considered domain. 



In all that follows, we consider stably stratified regions so that To^ - Tf > 0. 
We then non-dimensionalize the governing equations, taking our length scale to be the 



Stern 


1960 


Kato 


1966 



anticipated finger width, d = {kti^ / ga\Toz — Tf^\)^l'^ (Stern 
scale time, temperature, and composition as [t] = /kt, [T] = (Tq^ — Tf)d, and 
[p] = (a//3)(To2 — T'f)d. With the Prandtl number defined as Pr = v / and the diffusivity 
ratio as r = k^/ kt, the non-dimensionalized governing equations take the following form: 

V • u = 0, (6) 

^ + u ■ Vu^ = -Vp + (T - p)e, + V^M, (7) 

dT 

~di 
dp 



+ u • VT + w = V^T, and 



w 



rVV- (9) 
Equation ^ introduces a new parameter, Rq. This so-called "density ratio" is the 



ratio of the stabilizing entropy gradient to the destabilizing compositional gradient (Ulrich 
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T972|, 



V-Vad _ ajTo, - Tf) 



(10) 



If i?o < 1, the destabilizing compositional gradient dominates, and the system is unstable to 
the Ledoux criterion. One expects the system to be fully mixed by overturning convection. 
For Rq ^ 1, the stabilizing entropy gradient dominates, and the system is stable. The 



only possible transport is via diffusion. For intermediate values, 1 < i?o < 1/t, Baines & 



Gill (1969) showed that the system is unstable to fingering convection. It is this region 



of parameter space to which we now focus our attention. As in Traxler et al. (2011a), we 
define a reduced density ratio, 

r=5^, (IV 



that remaps the "fingering regime" to the interval of r G [0, 1] regardless of the value of r. 
We model fingering convection numerically, and solve Equations ^ through (|9| using a 



Traxler et al. 


(2011a|b 


) and 


Stellmach 



et al. (2011). To avoid spurious boundary effects, we apply triply-periodic boundary 



conditions on all perturbations T, fi, and u. This code uses no sub-grid scale model, so 
the scale on which energy is dissipated must always be fully resolved. In all simulations 
described below, unless specifically noted, the computational domain is taken to be a cube 
with side length lOOd. Since the wavelength of the fastest growing mode ranges from 6d to 
15d in the parameter regime considered in this paper, this implies that there are at least 5 
to 10 wavelengths of the fastest growing mode per lOOci. We test the required resolution 
by inspecting the compositional and vorticity fields, which are in general the least resolved 
ones given the small values of Pr and r selected. 



3. Results 



Traxler et al. 



(2011a 



,) explored a significant sample of parameter space, with Pr and 
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T ranging from 1/3 down to 1/30, as seen in Table [T| However, they only tested a few 
different density ratios in the cases with Pr, r ~ 1/30. Furthermore, owing to computational 
limitations, they did not explore the low Rq limit, the regime where we believe their 
prescription is most questionable (see § [T]). Here, we explore more comprehensively in 



Rq the parameters chosen by Traxler et al. (2011a) and also run simulations down to Pr 



r ~ 1/100. This enables us to test the validity of their prescription in more detail. It is 
computationally prohibitive to reduce Pr, r much lower than this, particularly at low Rq, 
where the simulations become increasingly turbulent]^ 



As in Traxler et al. (2011a), we report our results in the form of the Nusselt numbers. 



which are defined as the ratio of the total vertical flux to the diffused flux and expressed 
here in terms of the non-dimensional fields as El 



Nut = 1 - (wT), 
Nu„ = 1 - — (w/i). 

T 



(12) 
(13) 



The vertical turbulent fluxes, {wT) and (w/i), are measured by taking the average of the 



^We note that simulations at much lower values of Pr and r have been reported by 



Denissenkov (2010). Even though the latter are two-dimensional, it is unlikely that they are 



fully resolved. Whether under-resolved simulations yield satisfactory estimates for the fluxes 
is a different question that we shall address in a subsequent publication. The results from 



Denissenkov (2010) are consistent with ours within a factor of a few, which is expected when 



Stern et al. 



2001). 



comparing 2D and 3D results in flngering convection (see 

^The form in terms of dimensional variables w and T is Nu^ = 1 — {wT)/k,t{Toz — 
T^'^). Note that this definition of the Nusselt number describes the flux of the potential 
temperature, uriToz — T^'^), and not that of the absolute temperature. For a more complete 



discussion, see Mirouh et al. (2011). 
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products wT and wjj, over the entire domain. 



3.1. Typical and Atypical Simulations 

In Fig. [ij we present the temporal evolution of the thermal and compositional Nusselt 
numbers in a typical numerical simulation of fingering convection, taking Pr = 1/10, 
r = 1/30 and Rq = 3. This evolution begins with a period of exponential growth from t = 
to t = 160. The observed growth rate is found to be twice the growth rate of the fastest 
growing mode from the linear stability analysis of the governing equations, as expected. The 
transport peaks around t = 165, after which the fiux of both temperature and composition 
saturate (at a level slightly below the peak) and remain roughly constant for the remainder 
of the simulation. In order to visualize the saturation process, we present snapshots of the 
compositional field of the simulation at three stages before and after saturation in Figures 
[2]^a), (b), and (c), marked with vertical lines in Fig. [l] In Fig. [2|^a), which shows the state 
just as the primary fingering instability begins to develop, we see tall elevator modes with 
high n plumes sinking and low /i plumes rising. In Fig. |2](b), which is taken just before 
the peak in Fig. [T| we recognize the same vertical structures but notice that smaller-scale 
instabilities have distorted them. The latter destroy the elevator modes until the system 
achieves saturation in Fig. ^c). At this point, very little of the original elevator modes 
remain discernable, and quasi-steady, homogeneous turbulence dominates the dynamics. 

While most simulations behave exactly as Fig. [T| in some cases where Rq is chosen to 



be very close to the onset of overturning convection (i.e. for small r, see Equation (11)), the 
Nusselt numbers begin to increase gradually again after saturation. This behavior appears 
clearly in Fig. |3| for which Pr = 1/10, r = 1/30, and Rq = 1.1. We distinguish two phases 
post-saturation: from t = 100 to t = 600, the mean Nusselt numbers increase slowly and 
undergo quasi-periodic oscillations. We find that these oscillations occur with twice the 
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Fig. 1. — A typical example of the evolution of the thermal (top) and compositional (bottom) 



Nusselt numbers, Nut and Nu^ (see Equation (12)) in a simulation with Pr = 1/10, r = 1/30 
with Rq = 3. Time is measured in units of the thermal diffusion time across a length d. We 
note that the Nusselt curves grow exponentially, peak, and fall to a quasi-steady equilibrium 
state, which we call the saturated regime. The vertical lines indicate the times of the 
snapshots in Fig. [2] 
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(c) 

Fig. 2. — Snapshots of the compositional perturbation in the simulation shown in Fig. [T]at 
the characteristic times marked in Fig. [l] In Fig. [2]^a), prior to saturation, the composition 
field is dominated by tall "elevator mode" structures characteristic of linear fingering convec- 
tion. In Fig. ^h), though the original elevator modes are still recognizable, they have been 
disrupted by secondary instabilities. Fig. ^c) illustrates the "saturated regime" discussed 
in Fig. [T] Note that the elevator modes have been completely destroyed by the secondary 
instabilities, leaving a domain filled completely with homogeneous turbulence. 
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100 200 300 400 500 600 700 800 

t 



Fig. 3. — An "atypical" simulation showing the development and effect of large-scale in- 



stabihties on the Nusselt numbers, Nuy and Nu^ (see Equation (12)) in a simulation with 
Pr = 1/10, r = 1/30, and Rq = 1.1. We note that this simulation behaves quite differently 
from the one shown in Fig. [Tj After the initial saturation at t = 80, both Nusselt numbers 
gradually increase until t = 700. Both Nusselt numbers exhibit quasi-periodic oscillations 
at twice the buoyancy frequency (though such oscillations are more apparent in the thermal 
Nusselt number), which suggests the development of gravity waves. At t = 700, the trans- 
port increases dramatically in a manner similar to the onset of layers in the oceanographic 



case (e.g. simulations by Stellmach et al. 2011). The vertical dotted line marks the time 
at which we plot a snapshot of the compositional perturbation in Fig. |9l and the vertical 



dashed lines indicate the timesteps for which we calculate the density profiles in Fig. 10 
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buoyancy frequency, which is what we expect for gravity waves. Since hnear gravity waves 
do not lead to an increase in the mean flux, the slope of the mean Nut and Nu^ curves 
during that time indicate that the waves are highly nonlinear. At t = 600, the transport 
suddenly increases dramatically. Meanwhile, the mean density profile in part of the domain 



inverts, suggesting the formation of a fully convective layer (see § 5.2). This is reminiscent 
of the formation of gravity waves and layer development in oceanic fingering convection 



(Stellmach et al. 2011). We will discuss these large-scale instabilities and related transport 
properties in § [5j 



3.2. Extraction of the Fluxes 

To calculate the saturated fluxes in the homogenous phase of all simulations, i.e. prior 
to the onset of any large-scale dynamics discussed above, we use the following method. We 
identify the first local minimum in NuT(t) after saturation and set this to be the beginning 
of the saturated regime. We then fit the remainder of the time series to a linear function 
of time with a least squares fit, which yields a slope and a y-intercept. We use regression 
analysis to estimate an uncertainty for each parameter. If the uncertainty in the slope is 
greater than its magnitude, then we say that the Nusselt curve after saturation is fiat and 
set the end of the saturated regime at the end of the simulation. On the other hand, for 
a simulation like the one shown in Fig. [3| where the flux saturates initially but gravity 
waves begin enhancing the transport again shortly thereafter, we are only interested in the 
short interval just after saturation. If the linear fit is not sufficiently flat, we reduce the 
end time of the "homogeneous phase" progressively until the fitted slope is zero (within 
its uncertainty). Once the time interval for homogenous fingering convection has been 
identified, we calculate the saturated fiuxes by averaging them over this interval. The error 
is taken to be the root mean square of the fiuxes. 
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Table 1. Turbulent flux measurements and respective errors for low Pr, low r simulations 

of thermohaline convection 



Pr 


T 


Ro 


Nut 


(7t 


Nu^ 




Layers 


1/3 


1/3 


1.02003'' 


31.0 


4.0 


120.0 


10.0 


N 






1.08012'' 


23.0 


1.0 


92.0 


5.0 


N 






LIOOIS*" 


20.0 


1.0 


83.0 


5.0 


N 






1.2003^ 


13.5 


0.6 


60.0 


3.0 


N 






1.30045'' 


9.9 


0.5 


46.0 


2.0 


N 






1.50075'' 


5.7 


0.2 


28.3 


0.9 


N 






2.0015" 


2.16 


0.05 


9.5 


0.3 


N 






2.4021" 


1.3 


0.01 


3.5 


0.09 


N 






2.8027" 


1.037 


0.003 


1.33 


0.02 


N 






2.90285 


1.0159 


0.0009 


1.143 


r\ r\r\o 

0.008 


N 


1/3 


1/10 


1.01 


28.6 


0.4 


428.0 


6.0 


Y 






1.54" 


9.0 


0.2 


197.0 


4.0 


N 






3.97" 


1.62 


0.02 


38.7 


1.0 


N 






7.03" 


1.075 


0.003 


7.6 


0.2 


N 






9.46" 


1.0041 


0.0004 


1.41 


0.04 


N 


1/10 


1/3 


1.10015" 


8.0 


0.5 


33.0 


2.0 


N 






1.70105" 


2.16 


0.05 


8.7 


0.3 


N 






2.30195" 


1.24 


0.01 


2.9 


0.1 


N 






2.8027" 


1.019 


0.0007 


1.17 


0.006 


N 






2.90285" 


1.0065 


0.0003 


1.059 


0.003 


N 



1/10 1/10 1.003 11.36 0.09 169.0 2.0 Y 



-17- 



Table 1 — Continued 



Pr 


T 




Nur 


(7t 






Layers 






1.09^ 


8.3 


0.7 


143.0 


9.0 


N 






1.45^ 


4.3 


0.1 


86.0 


3.0 


N 






1.99^ 


2.5 


0.05 


54.0 


2.0 


N 






2.8^ 


1.73 


0.02 


34.8 


1.0 


N 






2.98'^ 


1.62 


0.01 


31.3 


0.6 


N 






3.25^ 


1.51 


0.02 


27.4 


0.8 


N 






4.96^ 


1.148 


0.003 


11.6 


0.2 


N 






7.03^^ 


1.032 


0.001 


3.9 


0.1 


N 






9.P 


1.0036 


0.0001 


1.36 


0.01 


N 


1/10 


1/30 


1.1 


5.6 


0.9 


540.0 


40.0 


Y 






1.5 


4.24 


0.06 


347.0 


5.0 


N 






2.0 


2.77 


0.02 


252.0 


2.0 


N 






3.0 


1.484 


0.008 


120.0 


1.0 


N 






6.0 


1.253 


0.004 


95.0 


1.0 


N 






8.0 


1.1382 


0.0006 


66.0 


0.3 


N 






10.963*^ 


1.059 


0.002 


32.4 


0.7 


N 






20.34^ 


1.0075 


0.0005 


7.1 


0.4 


N 


1/10 


1/100 


1.1b 


7.4812 




1584.517 




N 






1.3^ 


5.1979 




1285.765 




N 






1.5^ 


4.1962 




1124.055 




N 






1.7^ 


3.581 




1031.832 




N 
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3.3. Numerical Simulations 



A compilation of our new results with published data from Traxler et al. (2011a) and 



Radko & Smith (2012) are given in Table [T| Fig. |4] shows the corresponding values of the 
mean Nusselt numbers of both composition and temperature, in the homogenous saturated 



phase, as a function of the reduced density ratio r (see Equation (11)). As expected, we see 
that both Nut — 1 and Nu^ — 1 tend to as r — )■ 1. In contrast, as r — )■ 0, both Nu^ — 1 and 
Nu^ — 1 increase rapidly, presumably approaching values of Rayleigh-Benard convection in 
a triply-periodic domain ( [Calzavarini et aL 2005 Garaud et al. 2010). As an additional 
note, it is clear from Fig. Illthat Nu^ does not depend on Pr and r individually, but rather 



on their ratio. This effect had already been observed in Traxler et al. (2011a), and led them 
to propose their empirical scaling laws: 



Nut - 1 = Pr^/V^^/Vlr) 



Nu^- 1 




9[r) 



(14) 
(15) 



where g{r) and /(r) have the form ae '"'(I — rY, where for /(r), a = 264 ± 1, 
6 = 4.7 ± 0.2, c = 1.1 ± 0.1, and for g{r), a = 101 ± 1, 6 = 3.6 ± 0.3, c = 1.1 ± 0.1. 



As discussed in § [T| the prescription given by Traxler et al. (2011a) agrees well with 
the data for high Rq] however, these scalings underestimate the transport at low Rq, 
particularly as Pr and r decrease. This is illustrated in Fig. [5} in which we present the 



same data scaled by the suggested prescriptions from Traxler et al. (2011a). For the latter 



to remain an accurate representation of the data, all the points should lie on a single curve 
that depends only on r. It is clear from Fig. |5] that their prescription breaks down at very 
low r, particularly as Pr decreases. For this reason we now propose a theory that better 
captures the behavior at both large and small r. We derive this new theory with physical 
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Table 1 — Continued 
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^Data from Traxler et al. (2011a). 

'^Data from Radko & Smith (2012 ) with a domain size of 38.2 x 38.2 x 76.3 and a resolution 
of 256 X 256 x 512. The uncertainties of these fluxes are not know to us. 
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Fig. 4. — The Nusselt numbers as a function of the reduced density ratio, r, reported from 
the simulations in Table [l] (see references therein). Note that the compositional Nusselt 



number depends only on the ratio r/Pr and r, as noted in Traxler et al. (2011a). For all 
simulations, the error bars are smaller than the size of the plotted symbols. 
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Fig. 5. — This figure contains the same data presented in Fig. |4] scaled according to the 



model presented in Traxler et al. (2011a). Note that the scaling works well for large r, where 



the data collapse onto a single curve. However, at low r, the scaling breaks down. The 
points and line styles are the same as those in Fig. |4j 
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justification. 



Theoretical Model for Fluxes at Saturation 



Recently, Radko & Smith (2012) proposed a theory that characterizes the saturation 



of the fingering instability and predicts the thermal and compositional fiuxes at saturation. 
The fundamental idea of this theory is straightforward: the saturation is thought to be 
caused by a secondary instability (or "parasitic" instability) of the elevator mode and occurs 
when the growth rate of the secondary instability, a, is of the order of the growth rate of 
the elevator mode, A. The latter can be obtained by linearizing the governing equations 



and deducing the properties of the fastest growing mode. Radko & Smith (2012) then 



calculate the growth rate of the secondary instability using Floquet theory. Unfortunately, 
this procedure turns out to be computationally quite expensive. In stellar evolution codes, 
fiuxes need to be calculated for different values of Pr, r, and Rq corresponding to each 



radial mesh point considered at each timestep. Using the Radko & Smith (2012) method 
would be prohibitive for this application. 



We pursue a theoretical model that is effectively a simplified form of the Radko fc 



Smith (2012) theory but has the advantage of being much more straightforward and can 



thus be used in real time in stellar evolution codes. In what follows, we first review the 



derivation of the growth rate of the fastest growing mode. A, in § 4.1 (Baines & Gill 1969) 



In § |4.2[ we then propose a very simple analytical formula for the growth rate of the 
secondary instability, cr, and finally determine the Nusselt numbers expected when a ~ A. 
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4.1. Fastest Growing Modes 

The properties of the fastest growing fingering mode can be determined by linearizing 
the governing equations and assuming exponential solutions of the form 



for each of the velocity, temperature, composition and pressure perturbations (Baines & 



Gill 1969). This yields a cubic equation for the growth rate A, with coefficients that depend 
on the wave vector, k, and the governing non-dimensional parameters, Pr, r, and Rq. It 
can be shown that the fastest growing mode is independent of z, hence the term "elevator" 
mode. Since our equations are symmetric in x and y, A only depends on the magnitude / of 
the horizontal wavenumber and not on its direction. 



The resulting cubic then reads (Baines & Gill 1969): 



A^ + a2A^ + aiA + oq = 0, where (17) 
aa = /2(1 + Pr + r), 
ai = /^(rPr + Pr + r) + Pr ( 1 ^ 



ao = /VPr + /2pr ( r - ^ 



Ro) ■ 

To identify the fastest growing mode, we maximize A with respect to the wavenumber 
This yields a quadratic for A in terms of /: 

a2A^ + aiA + ao = 0, where (1^ 
as = 1 + Pr + r, 
ai = 2/2(rPr + r + Pr), 

ao = 3/VPr + Pr ( r ^ 
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Equations (17) and (18) can be solved numerically simultaneously to determine the 
wavelength and growth rate of the fastest growing mode. More interestingly for anyone 
interested in quick analytical estimates, it can be shown that in the limit of Pr, r <^ 1, 

^ if r < Pr < 1 

A^<; ■ , (19) 

^ if Pr < r < 1 



and 



A detailed derivation of these asymptotic limits (and higher-order terms) is presented in 
Appendix (Bl 



/2 ^ ^ (20) 
^1 + r/Pr 



4.2. Estimating the Nusselt Number from the fastest growing mode 



Radko & Smith (2012) consider all possible sources of secondary instabilities for the 
elevator modes. Here, for simplicity, we restrict our analysis to instabilities arising from 
the shear between adjacent elevators and neglect viscosity. As we demonstrate below, this 
yields very satisfactory results, and allows for a much simpler solution to the problem. 

The elevator modes are characterized by a vertically invariant flow fleld of the kind 

\J{x,t) = wqC'^^ sm{lx)ez = WE{t) s'm{lx)ez , (21) 

where A and I are the growth rate and horizontal wavenumber of the fastest growing mode, 
defined in the previous section. Without loss of generality, we have selected the horizontal 
phase of the mode to be sin(/x). This flow fleld is associated with a temperature and 
composition fleld 

Te(x, t) = Tqc'^* sin(/x) and /iE(a;, t) = /ioe'^* sin(/x) , (22) 
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where Wq, Tq and /io are related via (Radko & Smith 2012): 



ATo + Wo = -I'^To To 
\So + Rq^Wq = — r/^yUo =^ 1^0 = 



Wo 



A + /2 ' 

RqWo 



\ + tI 



2 ' 



(23) 
(24) 



since elevator modes are exact solutions of the full set of nonlinear equations (Equations 
(|6|-(|9])). At early times, the velocity shear, wsit), is weak and the elevator modes grow 



un-impeded. Using the definitions of the Nusselt numbers given in Equation (12) with 



Equations (21), (22), and (23), we find that 
NuT(t) = ' ^ 



LxLyLz J X Jy J z 




WE{t)TE{t) siv? {lx)dxdydz 



1 - -WE{t)TE{t) = 1 + 



Wn 



-,2\t 



2(A + /2) 

A similar expression can be obtained for Nu^(t). This shows that the Nusselt numbers 
initially grow exponentially with a growth rate that is twice that of the fastest growing 
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mode, as discussed in Section § 3.2, prior to saturation 



We now consider secondary shearing instabilities on the elevator mode fiow. Since 
Pr ^ 1, we assume that viscosity is negligible. In that case, the growth rate of the shearing 
modes, cr, is a simple increasing function of the velocity within the fingers. Since the latter 
increases exponentially, a also increases, while the growth rate of the fingering instability 
remains constant. We therefore expect that there will be a time where the two are of the 
same order of magnitude. At this point, the elevator modes are destroyed and saturation 
occurs. 



As in Radko & Smith (2012), we neglect the temporal variation of the primary elevator 



modes when evaluating the growth rate of the secondary shearing modes, and simply write 



U(a;, z) = We sin(Za;)e2 . 



(26) 
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It can then be shown, through dimensional analysis or through Floquet theory, that the 
shear growth rate is given by: 

a = KweI . (27) 

where is a universal constant of proportionality. By "universal," we imply that this 
constant is independent of any system parameter (r, Pr, Rq) or any property of the elevator 
mode. All information about the latter is contained in / and we- 

The dimensional argument goes as follows: the only relevant velocity in the system is 
that of fluid within the fingers, we, and the only relevant length scale is l/l, so the only 
relevant growth rate is weI- The same expression can be derived more rigorously using 
Floquet theory, as shown in Appendix [A} 

Assuming that saturation occurs when the shearing instability growth rate is of the 
order of the fingering growth rate can be written mathematically as 

a = KweI = cA, (28) 

where c is independent of the fundamental parameters of the system. This equation 
uniquely determines the velocity within the fingers at saturation to be 

= ^ (29) 

where C = c/K is another universal constant (to be determined). This expression is similar 



to the one used by Denissenkov (2010) to estimate the effective diffusivity by dimensional 
analysis. Our own model, however, departs from his analysis at this point and produces a 
result that more accurately fits the results of simulations. 

A given elevator mode with velocity we{x,z) = WEsm{lx) has a temperature profile 
Te(x, z) = TEsin(Zx) with Te = — xf^, for the same reasons as described above. Hence, at 
saturation, 
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using Equation (29). And by similar arguments, it can be shown that 



(31) 



where C is the same constant as in Equation (30) 



We compare this prescription with data from this study and those simulations from 



Traxler et al. (2011a) and Radko & Smith (2012) that approach the astrophysical regime. 



As illustrated in Fig. [6} for a value of C = 7 and r < Pr, we find that our theoretical results 
match all simulations remarkably well. It is important to note that, while C needs to be 
fitted, it is a universal constant and cannot depend on Pr, r, or Rq. The fact that we are 
able to use only one value of C to correctly fit the data suggests that the premises of this 
theory are correct. 

In Fig. [7], we show that — in the opposite case with r > Pr — we do not observe the 



same quality of fit. The full analysis by Radko & Smith (2012) is also unable to explain 



these results, so it is unlikely that that case is dominated by a parasitic instability like the 
one discussed in this section. Since stars are always in the former situation, the poorness of 
fit of the latter case is interesting but not relevant to our ultimate purpose. 



4.3. Asymptotic Expansions 



As such. Equations (30) and (31) yield the thermal and compositional Nusselt numbers 



provided A and I are known. Calculating these quantities requires the simultaneous solution 



of a cubic and a quadratic (see § 4.1), which can be done numerically quite easily. However, 
analytical approximations to Nut and Nu^ can also be obtained using asymptotic analysis 
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(see Appendix [B]) , to find that for r <^ (r, Pr) <^ 1, 



Nut - 1 ~ C'Pi (1 + 



+ 



Nu^-1 



+ 



1 + 



+ 



(32) 
(33) 



For the case with (r, Pr) <^ r ^ 1 and r/Pr 



Nut - 1 ~ C^Pr 



Nu^ - 1 



1 + 3 + 



r0 



1 + 



r 1 + 



+ 



(34) 
(35) 



Note that in this case, Nu^ depends only on and not on Pr or r individually, as was 



noticed by Traxler et al. (2011a). And finally for the case as r — )■ 1, (r, Pr) <^ 1, 



Nm 



1 ^ C'Fr 



2 i2 



27 



1 



9 ' 



Nu„ -l^C' 



Recall that the scalings for large r in 



:i-rY + 



27 



(1-0=^ + 



(36) 



Traxler et al. 



(2011a) were Pt^/^t^/^ and Pr^^V^^ 



/2 



which are close to the two latter cases. Note that the lowest order term for Nut — 1 
scales at best as Pr. In astrophysical objects, where Pr ^ 1, turbulent heat transport 
by homogeneous fingering convection is therefore negligible. Such is not the case for 



Nu,, 



1, which can become quite large, especially in the limit of r — )■ 0. We will discuss 



the implications of these results more in § |6} Finally, note that these asymptotic formulae 
should work well for the true astrophysical regime where Pr, r ^ 1. However, since none 
of our simulations are particularly close to this regime, the asymptotic approximation does 
not compare well with the numerical simulations we have run so far. 
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Fig. 6. — Comparison between the turbulent compositional flux observed in simulations 
(Table [T] and references therein), shown as symbols, to our theory (see eq. (31)), shown 
as curves, for the simulations with r < Pr. Since r is always less than or equal to Pr in 
astrophysics, the cases considered here are physically relevant. The value of the unknown 
constant C merely shifts the theoretical model predictions vertically by the same amount 
for all curves, and we flnd that the best flt has C ~ 7. Using this value, we flnd that the 
measured flux in simulations flts the model remarkably well. We rescale the fluxes by t/Rq 
to separate the theoretical curves for ease of viewing. 
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Fig. 7. — Same as Fig. |6]but for cases with r > Pr. The same quahty of fit is not observed 
in these simulations, and the difference between the data and the model is of the order of a 
few. By contrast with Fig. 16} We rescale the fluxes by r / Rq to be consistent with Fig. [6] 
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5. The Development of Layers 

As discussed in § [T| substantial evidence from oceanography exists to suggest that 
mixing rates beyond the levels discussed in the previous section are possible. In fingering 
regions in the ocean, thermohaline staircases — layers of overturning convection separated by 
thin fingering interfaces — have been observed and have been linked with dramatic increases 



in the transport of composition and heat (e.g. Schmitt et al. 2005). Thus, it is important 
to consider whether layers can form in the astrophysical case, and whether they also lead 
to increased transport. 



5.1. The 7-instability 



Radko (2003) developed the so-called "7-instability" theory to explain the formation 



of thermohaline staircases in fingering regions of the ocean. The theory was found to 
agree with direct numerical simulations of staircase formation in two and three dimensions 



(Radko 2003; Stellmach et al. 2011) and has since become well-accepted in the physical 



oceanography community. 



The original 7-instability theory is a linear mean-field instability theory. Radko 



(2003) first averaged the governing equations of fingering convection (see Equations ^ 
through ([9])) spatially, over length scales that span several fingers, and temporally, over 
multiple finger lifetimes. He then performed a linear stability analysis of the averaged, 
or "mean-field," equations. He argued that the stability of the system depends on the 
quantity 7turb, defined as the ratio of the turbulent thermal fiux, Nut — 1, to the turbulent 
compositional fiux, r(Nu^ — l)/Ro. More specifically, he showed that homogeneous fingering 



convection is unstable to layering if and only if d'jturh/dRo < 0. Traxler et al. (2011a) 
later improved this theory and showed that in the astrophysical case, this result still holds 
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as long as 7 is redefined as the ratio of tlie total fluxes, diffusive and turbulent]^ in other 
words, a system is unstable to layering if and only if 

^^<0with7 = ^. (37) 



dRo rNu 



As mentioned in § 3^, some of our simulations at very low r exhibit properties 
consistent with layering. In order to determine whether these layers form via the 
7-inst ability, we plot in Fig. |8] the variation of 7 with Rq for various values of Pr and r. 
Unlike 7turb, which does decrease for some r at low Pr, r, the quantity 7 always seems to 



increase. As discussed by Traxler et al. (2011a), this is because Nut and Nu^ are both 



dominated by diffusive fluxes at low Pr, r and the corresponding "diffusive-7," Rq/t, 
increases with Rq. In short, since 7 is always an increasing function of Rq in our simulations, 
the observed layers cannot be forming by the 7-instability. Thus, while the 7-instability 
can lead to the development of staircases in high Pr oceanic fingering, it is unlikely to be 



relevant for astrophysical fingering, confirming the results of Traxler et al. (2011a) 



5.2. Simulations with Low Rq 

Despite the fact that 7 is a strictly increasing function of r, we do observe layer 
formation in some simulations. We now study these results in more detail, focusing on 
the simulation presented in Fig. [3| which has Pr = 1/10, r = 1/30, and Rq = 1.1, and 
shows strong evidence for the formation of a convective layer. Fig. |9] shows a snapshot of 
the compositional perturbation at t = 825; large-scale plumes characteristic of overturning 
convection are clearly visible. We can check to see whether these plumes are associated 



^The discussion of Denissenkov (2010) on layer formation uses 7turb rather than 7, which 
Traxler et al. (2011a) showed to be incorrect for the astrophysical case. 
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Fig. 8. — The total flux ratio, 7, using the same data as in Fig. |4} Note how 7 increases 
monotonically with Rq, which implies that thermocompositional staircases cannot form by 
the 7-instability. 
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Fig. 9. — A snapshot of the compositional perturbation from the same simulation as in 
Fig. |3] (Pr = 1/10, r = 1/30, and Rq = 1.1) at t = 825. Large convective plumes are clearly 
visible and span a significant fraction of the compositional domain. The interface is harder 
to identify but can be seen near the bottom of the plot, where finger structures are still 
apparent. Note that it is far from "flat," and instead is highly distorted by the large-scale 
plumes in the convective layer. 
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Fig. 10. — The horizontally-averaged density profile, (ptot) (see eq. (38)) as a function of 
z for the simulation shown in Fig. |3] at the four times marked by vertical dashed lines. 
The background density is shown for reference as a solid line and is always decreasing with 
height. The sharp density transition in the lower part of the domain is the interface. At 
times t = 725, 750, and 800, there exist regions over which the density increases with height, 
which we call a density inversion. The only time without an inversion, t = 775, corresponds 
to the point where the Nusself number is at its minimum (see Fig. [3]). 
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with an inversion of the horizontally-averaged total density profile, which is given by 

{p,^,) = -(l-R-^)z-{T) + {fi). (38) 



Fig. 10 shows (ptot) as a function of height for the simulation shown in Fig. [3] at four different 
times selected during the later stages, where the Nusselt numbers increase significantly. 
The solid line indicates the linear background gradient for reference. Note how each of 
the four density profiles shows a sharp, stably stratified transition region between z = 25 
and z = 45, the interface. Three of the four profiles also have a region where the total 
density increases with z. This demonstrates that the homogeneous fingering has indeed 
transitioned into a layered system with convective layers separated by stable interfaces. At 



t = 775, shown in Fig. 10, the convective layer is just beginning to develop, and a slight 
inversion of the density profile appears in Fig. [Toj As the simulation approaches the stage 
of maximum transport at t = 750 as seen in Fig. [3} the density inverts more strongly. The 
ensuing period of active mixing thickens the interface again and flattens the density profile. 
At t = 775, the convective mixing briefly shuts down. This process repeats, as can be seen 
at t = 800, when the density profile once is beginning to invert. Note that the interface 
appears to move up and down during this period, which we attribute to the advection of a 
large quantity of compositionally dense material by convective plumes. 

This is not the only such simulation in which we observe the formation of staircases: 
we see this behavior at low r (< 0.003) for several values of Pr and r, including Pr = 1/10, 
r = 1/10 and Pr = 1/3, r = 1/10. In all these cases, the Nut and Nu^ time series as well 
as the properties of the layered state are qualitatively similar to the one shown in Fig. [3} 
Fig. m and Fig. [10 



Having established that this layering transition cannot be due to the 7-instability, we 
look instead at the original mechanism proposed by Stern ( 1969[ ) (see also [Stern et al.||200lD 



for the formation of thermohaline staircases in the ocean. In addition to the 7-instability, 
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fingering convection is also known to excite large-scale gravity waves through another 



mean-field instability, the "collective instability." Stern (1969) studied this effect, and 
showed that large-scale gravity waves can be excited under some circumstances and grow in 
amplitude exponentially. These waves transport momentum, heat, and heavy elements until 
they begin to "break," at which point the advected material mixes with the background. 
He went on to suggest that if the waves break on large enough scales, this could cause 
enough local mixing to trigger the formation of staircases. 

In direct numerical simulations that have been run at oceanographically relevant 
parameters (Pr = 7, r = 0.3), [Stellmach et al. (2011) found that gravity waves are indeed 



excited by the collective instability but do not trigger layers. Layer formation in their 
simulations — and thus presumably in the ocean — is caused by the 7-instability rather than 
the collective instability. Our simulations suggest, however, that the converse may be true 
in the astrophysical parameter regime. As discussed in Fig. |3| gravity waves are observed 
to dominate the dynamics of the system between t = 100 and t = 600. Furthermore, the 
gradual increase in Nut and Nu^ during that phase suggests that the waves have high 
enough amplitudes to interact nonlinearly with each other, and with the background, by 
contrast with the oceanographic case. Indeed, linear waves cause oscillations in Nuy and 
Nu^ without changing their mean values — only nonlinear waves can affect the mean. Since 
the waves in Fig. |3] are clearly strong enough to interact nonlinearly, we can hope that 
they may also break and cause mixing on larger scales, following the original hypothesis 



proposed by Stern (1969). 



Checking this hypothesis in detail is difficult with the available simulations. However, 
we can at least determine whether our observed waves do indeed appear in accordance with 
the collective instability theory. A system is unstable to the collective instability if the 
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Fig. 11. — Stern number as a function of r, calculated from the data presented in Fig. |4] 
(and Table [l| using equation eq. ([39|)). Those simulations that develop into layers are denoted 
with a large black circle. A Stern number exceeding order unity suggests the possibility of 
triggering of gravity wave growth but doesn't necessarily imply layer formation. Here we 
see that the only cases in which layers form occur when A > 10^, which is marked with a 
horizontal dotted line. In addition, not all such simulations develop into layers. This could 
be due to the domain size or integration time of the simulation or perhaps implies that the 
condition for layer formation also depends on r. 
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Stem number, in our notation defined as 

(Nut - 1) (Ttuib " l) 



A 



Pr (1 - i?o ^) 



(39) 



exceeds order unity (Stern et al. 2001). Note that this expression uses 7turb 



i?o(NuT — l)/r(Nu^ — 1) instead of 7. For the simulation shown in Fig. and Fig. 10 



Nut — 1 = 6.55 and 7turb = 0.438, so that the Stern number is 924. This implies that the 
system is indeed unstable to the collective instabihty. 

More generally, we find that all the simulations in which we observe layer formation 
have Stern numbers of the order of 10'^ or more. However, not all such simulations go into 
layers, as can be seen in Fig. [TT| so it is likely that the condition for layer formation through 
the collective instability also depends on the value of r. It is also possible that these gravity 
waves grow in all simulations with A > 1 but that the growth rate depends on the Stern 
number, resulting in undetectable growth if the Stern number is too low. If this is the case, 
it will be difficult to see layer formation in models with smaller A because these simulations 



are expensive and difficult to integrate for long. We discuss this more in § 6.3 



Much work still needs to be done to identify the parameters for which we expect layer 
formation in the astrophysical regime and to quantify transport in the layered case (see 
Wood et al. (2013), in prep, for related efforts in the case of semi-convection). We have 
only three simulations that develop into layers, and we require a more extensive sample 
of simulations at low r in order to better understand the conditions of layer formation 
and their properties as a function of Pr, r, and r. In every simulation that we have run 
that exhibits layers, a single convective plume spans the majority of the domain, and 
we do not see the development of multiple layers, as has been seen in simulations of the 



fingering regime in the oceanographic case (e.g. Stellmach et al. 2011) and in the case of 



semi-convection (Rosenblum et al. 2011 Mirouh et al. 2012). To characterize this process 



we must increase the size of the domain and extend the integration time to let the layers 
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develop more naturally. This will allow us to describe more completely the size scales 
and transport of thermocompositional layers in astrophysical fingering convection. This 
problem, however, is beyond the scope of this paper and is deferred to a later publication. 



6. Discussion 

In the preceding sections, we have used the results from § [3] to formulate a theory 
that models the Nusselt numbers of homogeneous fingering convection in the astrophysical 
regime. We now provide guidance on how to convert real stellar parameters to the 
non-dimensional parameters for use in one-dimensional stellar codes. We quantify the 
range of parameters where layer formation via the collective-instability may be possible. To 
provide perspective on the impact of these results, we conclude by discussing the prediction 
from this theory in several astrophysical systems. 



6.1. A Method for Applying the Model 

The first step toward using our model in stellar evolution calculations is to identify the 
non-dimensional parameters of the system in question at each grid point and/or timestep. 
Recall that Vi = v / k^t and r = n^/ k.^- The expressions for the microscopic viscosity and 



diffusivities can be found in many textbooks of plasma physics (e.g. Chapman & Cowling 



1970). In stars, the value of Pr typically ranges from 10 ^ to 10 ^ and that of r from 10 



down to 10 ^, with Pr > r. 



The density ratio, Rq, is given by 



where V = (dlnT/dlnP) is the actual temperature gradient, Vad = {dlnT/dlnP)s is the 
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adiabatic temperature gradient, = [din fi/dlnP) is the mean molecular weight gradient, 



and 6 = —{d\np/d\nT)pf^ and = {d\n p/d\n fi)pT (see Kippenhahn &; Weigert 1990, for 
a detailed description of these gradients). The quantities V and can be calculated from 
the stellar model grid. The other values, Vad, 0, and 6, can be calculated from the equation 
of state. Note that Rq can take a broad range of values in stellar interiors, but the fluid is 
fingering-unstable only when 1 < i?o < t^^- 

Once Pr, r, and Rq are known, one can then calculate Nut and Nu„ according to 



Equations (30) and (31) for a numerically precise answer or by the asymptotic expression 
in § 4.3 if a rough estimate is all that is needed. Since Nu^ is the total flux of composition 
divided by the background diffusive flux, the effective compositional mixing coefficient for 
modeling the transport of fingering convection is given by 

= Nu^K^. (41) 

An analogous expression can be derived for the heat transport but this term is nearly always 



negligible (see § 4.2) when Pr ^ 1 unless layers form. We ignore it from here onward. By 
contrast, the turbulent compositional transport can be significant, as can be seen in Fig. [12 
where we plot the logarithm of Nu^ — 1 as a function of r/Pr = k^/z/ and r, assuming 
Pr = 10^^. Near the marginal stability limit, r = 1, the compositional turbulent transport 
is negligible as well, and Nu^ — 1 drops to zero. However, for low r, and particularly for 
low r/Pr, the turbulent compositional transport increases by many orders of magnitude. 
This plot is intended as a quick estimate for a large range of parameters; more accurate 



estimates can be obtained from the expressions given § 4.2 in § 4.3 



The results presented so far concern only the case of homogeneous fingering convection. 



As discussed in § |5.2[ however, layer formation seems to be possible in some cases and leads 
to significant increases in transport about the homogeneous levels. To determine when 
this is likely, we use the theory from § |4] to estimate the Stern number as a function of 
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Fig. 12. — Color plot illustrating the non-dimensional turbulent compositional transport as 
a function of k^/i^ and r from the model presented in § 4.2, We assume that v / = 10^^ 
for illustrative purposes. This plot serves as a quick reference for those wanting an order of 
magnitude estimate of the chemical transport for given k^, z/, and r. Note that the transport 
becomes dominated by diffusion near r = 1, as expected. For low r, particularly at low k^/ z/, 
the transport increases by five to six orders of magnitude. 
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Fig. 13. — Color plot illustrating the Stern number as a function of i^^/v and r for the same 
data as in Fig. 12 The green contour marks the 10^ line. If the parameter regime of layers 
is limited to A> 10^, the only viable range for layer formation is r < 10~^; however, it may 



be possible for layers to form in other cases, see discussion in § 6.3 
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"/Pr = ti^/v and r and show the results in Fig. 13 The Stern number is largest for small r 



and large r/Pr. Recall from § 5.2 that we observe layer formation for a Stern number near 
or above 10'^. Using the same criterion here, which is the best that can be done without 
further information, suggests naively that layer formation will only be relevant in stellar 
interiors (where r/Pr < 0.1) for r < 10~^, which is so close to actual overturning convection 
that it may in fact be irrelevant for practical purposes. Of course, as mentioned in § |5.2 , 



much work still needs to be done to confirm these results. This will be the subject of a 
subsequent pubhcation. 



6.2. Impact of the New Model on Astrophysical Applications 



We now reexamine several problems that have been recently discussed in the context of 



fingering convection, namely planet pollution (Vauclair 2004; Garaud 2011) and the impact 
of ^He fusion in red giant branch stars (Charbonnel & Zahn 2007 Denissenkov 2010). 



As first discussed by Vauclair (2004), fingering convection determines the timescale 



during which planetary material remains detectable on the surface of a star after infall. 



Garaud (2011) used the fingering transport coefficient from Traxler et al. (2011a) to 



determine that this timescale is too short for the excess metallicity of stars hosting 



planets to be related to planet infall. Since the prescription by Traxler et al. (2011a) only 



underestimates the transport at low r while still recovering the scalings at high r, the use 
of our new, improved model will only increase the turbulent mixing rates, which serves to 



decrease the mixing time, so the qualitative conclusion of Garaud (2011) still holds: the fact 



that planets are more readily found with host stars of high metallicity must be a primordial 
effect. 



To address the case of red giant branch stars, we compare our results to the numerical 
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simulations of Denissenkov (2010). For his choice of governing parameters, Pr = 4 x 10~^, 
r = 2 X 10^^, and Rq = 1700, he finds Nu^ = 998. He concludes from this result that 
fingering convection cannot provide all the additional mixing needed in red giant stars. Our 
model finds a remarkably similar value of Nu^ = 1294 for the same parameters, and thus 
we confirm his conclusion that additional mixing is necessary to explain observations. 



6.3. Summary and Future Work 



We have developed a theoretical model for the vertical transport of heat and 
composition by homogeneous fingering convection that reproduces the results of numerical 
simulations from this study and from previous studies. We have found that the turbulent 
heat transport is negligible, but that compositional transport can be important. In 
particular, this model predicts more transport for weakly stratified systems than was 



initially suggested by Traxler et al. (2011a). This model can be extrapolated to stellar 



parameter regimes and used as a simple method of calculating the compositional turbulent 
diffusion coefficient in fingering regions of stellar interiors and can be done in real time in a 
stellar evolution code. 

We have also found the first evidence for layer formation by fingering convection in 
the astrophysical parameter regime. We have identified the collective instability as the 
origin of layer formation. By contrast, in the oceanographic and semi-convective cases, the 



7-instability causes layer formation (e.g. Stellmach et al. 2011 Mirouh et al. 2012). Much 
work still needs to be done to characterize the conditions for layer formation and transport 
by layered convection in this case. With our limited computational resources, we have found 



that layer formation appears to require high Stern numbers (see Equation (39)) to occur. If 



true, this would imply that layers are only likely to form in regions that are very close to 
being fully convective anyway. However, it remains to be seen whether lower Stern numbers 
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can also result in layer formation. To study this will require additional simulations and 
longer integration times than have been covered in this paper. In addition to the conditions 
of layer formation, it is also important to develop a theoretical or empirical model for the 
transport of this case. Our simulations show that transport increases significantly when 
layers form. Concurrent work by Wood, Garaud, and Stellmach (2013) reveals that thermal 
and compositional transport in the layered semi-convective case scales with the layer height 
to the power of 4/3. If this scaling also applies here, we may expect that very significant 
transport rates could occur for large enough layer heights. What determines the latter, 
however, will also require new theories. 
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A. Floquet theory for the secondary shearing instability of fingers 

As discussed in § |4| dimensional analysis can be used to show that the shearing growth 
rate a of the fingering elevator modes must be proportional to the vertical velocity within 
the finger, times their horizontal wavenumber, I. The same conclusion can be reached 
using a more formal approach through Floquet theory, detailed in this Appendix. 



We loosely follow the method described in Radko & Smith (2012), but use a number of 



additional simplifications that enable us to obtain analytical solutions. 

We first restrict our analysis to 2D fiows. This simplifies the problem greatly, and it 
can be shown that the final result (i.e. a oc wl) would be the same in 3D. Secondly, we only 
consider the effect of shear in driving secondary instabilities, and neglect that of buoyancy 
(i.e. we neglect temperature and salinity perturbations). The 2D elevator modes are then 
well-described by the following velocity field: 



U(x,t) = w-^it) sin(/x)e2 . 



(Al) 



The choice of the phase of U(x,t) (i.e. sine or cosine) is arbitrary. We choose the sine for 



consistency with Radko & Smith (2012). As in Radko & Smith (2012), we then ignore the 



time- dependence of the elevator mode velocity field U, setting w-E,{t) = we where we is 
constant. Finally, we neglect the effect of viscosity entirely (which is justified in the low 
Prandtl number astrophysical limit). 

We then consider perturbations u' = {u',0,w') to this basic fiow, driven by secondary 



shearing instabilities. As in Radko & Smith (2012), we work with a stream function ip 
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defined so that u' = —dip/dz and w' = dtp/dx. The hnearized momentum equation reads: 

du' 

— + u • VU + U ■ Vu = -PrVp, (A2) 
at 

where the nonlinear term U • VU naturally vanishes, since the elevator modes are fully 
nonlinear solutions of the governing equations. Written in terms of the stream function, 
this becomes: 

— ( V V) + WE sm{lx) — ( V^^) + weP sm{lx) ^ = 0. (A3) 

This is a linear partial differential equation for ip, with coefficients that are constant in 
time and in the coordinate z. We may write the solutions as 

V^(x,2,t) =^(x)e*'""+"* (A4) 

to satisfy periodic boundary conditions in z, as in the original problem. Floquet theory 
further shows that the x-dependence of the solution must be in the form 

N 

^ = e^fi- J2 V^ne*"'", (A5) 

n=-N 



where we have used for the sake of clarity the same notation as Radko & Smith (2012) 
for the Floquet term e*-^'^, where / is the Floquet coefficient. Since / needs not be an 
integer, the shearing modes are usually quasi-periodic rather than periodic. Note also that 
the summation term should in reality be taken from n = —oo to +oo, but is written here 
already in its approximate truncated form, in view of a numerical implementation of the 
problem. 



Plugging this expression into Equation (A3), and projecting the result onto individual 
Fourier modes, we obtain the following tri-diagonal linear system: for each value of n 
ranging from — to A^, we have 

[m' -f + l\f + n - If] - a^^ [m' + l\f + nf] + 

^^„H_i - + nf + n + If] = 0, (A6) 



with the imphcit assumption that ^/'-at-i = ^n+i = 0. Each secondary mode of instabihty 
is identified by the pair {f,m), and grows with rate a{f,m) obtained by setting the 



determinant of the hnear system (Equation (A6)) to zero. As in Radko & Smith (2012), 
we can then scan over all possible values of m and / to find the growth rate of the 
fastest-growing secondary instability. 

The beauty of this method is that, by contrast with the more general system studied 



by Radko & Smith (2012), the scalings for a can be obtained without solving for it at 
all! Let's define the rescaled vertical wavenumber m = m/l and the rescaled growth rate 
a = a /w-eI. The system of equations described by Equation ( |A6[ ) for n = —N..N can then 
be re-cast in the form 

-li^n-i [rh'-l + if + n- ir] - [rh' + (/ + n)'] + 

^^„+i[m2-l + (/ + n + l)2] = 0. (A7) 

The rescaled growth rate a is now a function of m and / only, so that maximizing a over all 
possible values of m and / yields one universal constant. In other words, the fastest growing 
secondary shearing mode satisfies a = K, where i^' is a universal constant, independent of I 
or We- This finally implies, as discussed in the main text, that 

a = KweI, (A8) 

where the value of K is somewhat irrelevant since it can be folded into the constant C 
defined in § |4j 



B. Asymptotic Analysis 



In Section § |4| we derived an analytical expression for the thermal and compositional 
Nusselt numbers (see (30) and (31)), in terms of the growth rate A and wavenumber / of the 
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fastest growing elevator mode. The latter can be found semi-analytically by solving a cubic 



and quadratic simultaneously, given by Equation (17) and (18) 



While this can be done exactly with a Newton-Raphson relaxation method (or any 
other numerical method), one may be interested in a quick fully analytical estimate of 
these fluxes instead. In this Appendix, we derive an approximate formula for the turbulent 
fluxes, which is increasingly accurate in the asymptotic limit of Pr, r — )■ 0. To do so, we 



perform an asymptotic analysis of Equation (17) and Equation (18) in these limits. Since r 



is always of the order the Prandtl number in most astrophysical applications, we simplify 
our analysis by considering the two limits at the same time, defining the quantity = r/Pr, 
and requiring it to be of order unity. 

Asymptotic analyses are always much easier to do if one has a good idea of the behavior 



of the exact solution. The numerical solutions to Equation (17) and Equation (18) are 
shown in Fig. [M] (see data points). We see that the growth rate of the fastest growing mode 
A as a function of the reduced density ratio r follows different asymptotic laws depending 
on the range of r considered. For r ^ (Pr, r), A appears to be constant and proportional 
to \/Pt. For r <^ 1 but not in the previous limit, A appears to be proportional to Pr, and 
decreases as r^^/^. Finally, in the limit of r — )■ 1, A drops to 0. We now investigate these 
three limits in more detail. 



B.l. First Regime: r ^ (Pr, r) ^ 1 

In the limit of r ^ (Pr, r), we find numerically that A is independent of r. The same 
is true for the wavenumber I of the fastest growing mode. This suggests that both A and 
I are continuous in the limit r — 0, and that one may study this first regime simply by 



solving Equation (17) and Equation (18) with r = (alternatively, Rq = 1). The resulting 



equations are (using the definition of (j) given above) : 



A' + + Pr(l + <i))) + ArPr(l + (/)(Pr + 1)) + m> + m(Pr(/) - 1) = , (Bl) 

and 

(1 + Pr(l + + 2A/2pr(l + 0(Pr + 1)) + S/^Pr^c/) + Pr(Pr0 - 1) = . (B2) 

Subtracting P times tlie second equation from the first yields the simphfied cubic: 



A^ - ArPr(l + 0(Pr + 1)) - 2rPr> = 



(B3) 



Next, we study the behavior of A and / with varying Pr and r (via cj)). We observe that 
the exact solution for A varies as Pr^/^ while P is more-or-less independent of Pr. Using 
this information, we propose the following asymptotic expansions: 



A = Pri/2 (^Ao + AiPri/2 j + . . . ^nd = /2 + /2p^i/2 ^ . . . 



(B4) 



expanding P rather than / since Equation (17) and Equation (18) only depend on the 
former. 



Plugging these two ansatz into Equation (B2) and (B3), and solving for all unknowns 
at their respective orders in Pr yields 



A : 

/2 



Pr - VPr^ + rPr + ■ 
1 



Vp^ 1 



rPr 



(r + Pr)2 



+ ■ ■ ■ 



(B5) 



As can be seen in Fig. 14, this approximation does well for small Pr and r for 



r < (Pr,r). Note that although the only illustrated cases are those with Pr = r, we have 



tested cases down to 



10-2. 



We may then use Equations (30) and (31) to estimate the thermal and compositional 



- 52 - 



Nusselt numbers. Keeping only the leading order terms yields 



Nut - 1 ~ [Ft (1 + 0) - Pr^^^ 



3 + 40 + 



Nu^- 1 



2 + 20 + ^^2 



+ 



Pr 



1 + 



+ 



(B6) 
(B7) 



This too is compared to the numerical solution in Fig. 15 



B.2. Second regime: (Pr, r) ^ r <^ 1 



In this regime, since (Pr, r) ^ r, we expand Equation (17) and Equation (18) first in 
Pr and then in r. Inspection of the numerical solutions suggests that A is proportional to 



Pr, while is more-or-less independent of it. Seeking solutions of Equation (17) of the form 
A = PrA, to the lowest order in Pr, yields 



A[/^(0+l)-l] +2/V = 



(B8) 



which can be easily solved to find 



2/6 



/4(l + 0) - 1 



(B9) 



Plugging this into Equation (18), keeping only the lowest order terms in Pr, and collecting 



the results in terms of I yields 
-/i2(l-0)' + /«(0+l) 



+ 1 



+ 1' 



1-20 +-(0 + 1) 



r — 1 



(BIO) 



which is a cubic equation for Unfortunately, this cubic has no simple solution, so we look 
for an expansion in the other asymptotic variable, r. Inspection of the numerical solutions 
suggests that the expansion needs to be in terms of r^^'^, so we set = /q + Ify/r + l2r H . 



Using this ansatz in Equation (BIO) yields, to lowest order 



/«(l + 0)2-2/o^(l + 0) + l = O, 



(Bll) 



- 53 - 



which we can solve to find = 1/(1 + This, on its own, causes A to diverge (see 



Equation (B9)), so we need to go to the next order in ^/r. We then find that 



1 



- 2- 



r0 



1 + 0)5/2 



+ 



which eventually yields 



A 



Prr / Pr 

3r 



3/2 



+ 



(B12) 



(B13) 



r \ Pr + r ^ 

keeping the first two terms in the expansion only. As expected from the numerical solutions 
discussed above, A scales like r^^/^Pr to lowest order. 



As in the previous section, we use Equation (30) and Equation (31) along with these 



expansions to find approximate expressions for Nuy and Nu^: 



Nut - 1 ^ C^Pt\- ' ^ + 

r 



-2, 



r (1 + 



+ 



Nu^ - 1 ^ C 



1 + 3 + 



r0 1 + 



+ 



(B14) 
(B15) 



We compare the results of the asymptotic expansion to the numerical solution in 



Fig. 14 and Fig. 15, where we look at A and Nu^ — 1, respectively. We choose to set the 
transition between the two regimes at r = (r, Pr). For those simulations where r is several 
orders of magnitude smaller than Pr, we recommend choosing this point at r = r. It is 
clear that the asymptotic expansion does an excellent job of reproducing the true numerical 
solution in the limit considered. We also see that, as suggested by the numerical data, Nu^ 



only depends on rather than on Pr or r individually as was noticed by Traxler et al. 
(|2011aD. 
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Fig. 14. — Comparison of the asymptotic expansions in Equations (B5) and (B13) to the 
numerical solution of A for small r using C = 7 keeping up to the second order. The 
asymptotic expansion switches between two regimes of low r, r ^ (r, Pr) and r ^ (r, Pr), at 
r = (r, Pr). The data points represent the numerical solution and the curves represent the 
asymptotic expansion. Note that the aysmptotic expansion fits the data remarkably well in 
the respective limits considered. 
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Fig. 15. — Comparison of the asymptotic expansions in Equations (B6) and (B14) to the 
numerical solution of Nu^ — 1 from Equation (31) for small r using C = 7 keeping up to the 
second order. Note that this retains the same quality of fit as Fig. 14 
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B.3. Third regime: r — ?■ 1 



Equations (B9) and (BIO) were derived witliout any assumptions on the value of r, 



aside from the fact that it needs to be much larger than Pr and r. We can therefore 
use them directly to study the limit of r — )• 1. To do so, we set e = 1 — r, and rewrite 



Equation (BIO) in terms of the small parameter e: 



P(i 



+ + 1) [2 + (0 + l)e)] - [3 + 2e(0 + 1)] + e(l + e) = 0. 



(B16) 



Inspection of the solutions for l'^ in the limit of small e suggests the expansion 



/■^ = e(/g + elf + ■ ■ ■). Plugging this into Equation (B16) then yields, to lowest order 



I' 



1 



(BIT) 



Plugging this expression into Equation (|B9|), we get 



A = 2r 



-1 3/2 



-(1-r) 



'1 — r' 



(B18) 



As before, we use Equations (B17) and (B18) with Equations (30) and (31) to find 



expansions in Pr and r for Nut and Nu^ in the limits Pr — i- and r — 1. We find 

'4 



Nut-1 ^CW(j)'^ ( -d-r) 



27 



1 + 



1-r 



Nu,, 



1 ~ d - (1 



27' 



:i-r) 



(B19) 



We compare the asymptotic expressions to the numerical solutions for A and Nu^ in 



Fig. 16 and in Fig. 17 The expansion clearly works best for r > 0.5. Note that these 



asymptotic expansions are designed to work only when ~ 1 and should be used with 
caution if 6 <^ r. 
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Fig. 16. — Comparison of the asymptotic expansion in Equation (B18) to the numerical 
solution of A as r —7- 1. The data points represent the numerical solution and the curves 
represent the asymptotic expansion. Note that the aysmptotic expansion fits the data well 
for its intended regime. 
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Fig. 17. — Comparison of the asymptotic expansion in Equation (B19) to the numerical 
solution of Nu^ — 1 from Equation (31) as r — t- 1 with C = 7. Note that this retains the 
same quality of fit as Fig. [16 
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